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Abstract 

Spatial adaptation procedures for the accurate and effi- 
cient solution of steady and unsteady inviscid flow problems 
are described. The adaptation procedures were developed 
and implemented within a three-dimensional, unstructured- 
grid, upwind-type Euler code. These procedures involve 
mesh enrichment and mesh coarsening to either add points 
in high gradient regions of the flow or remove points where 
they are not needed, respectively, to produce solutions of 
high spatial accuracy at minimal computational cost. The 
paper gives a detailed description of the enrichment and 
coarsening procedures and presents comparisons with ex- 
perimental data for an ONERA M6 wing and an exact so- 
lution for a shock-tube problem to provide an assessment 
of the accuracy and efficiency of the capability. Steady and 
unsteady results, obtained using spatial adaptation proce- 
dures, are shown to be of high spatial accuracy, primarily in 
that discontinuities such as shock waves are captured very 
sharply. 


Introduction 

Considerable progress in developing methods of dy- 
namically adapting computational meshes based on the nu- 
merical solution of partial differential equations has been 
made over the past decade. 1 These methods are being de- 
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veloped to produce higher spatial accuracy in such solutions 
more efficiently. Spatial accuracy is obviously important 
when modeling continuous equations with a discrete set of 
points. It is generally understood that accuracy is improved 
when the number of mesh points in a fixed computational 
domain is increased. Associated with an increase in the 
number of mesh points, however, are increased computer 
run times and memory costs. Hence, for efficiency, it is 
important to enrich meshes locally based on the numerical 
solution, in contrast to using globally fine meshes, to min- 
imize the total number of mesh points and hence minimize 
the cost for a given spatial accuracy. The methods of mesh 
refinement can be separated into three general categories: 
(1) mesh regeneration, (2) mesh movement, and (3) mesh 
enrichment. 

The first method, mesh regeneration, places the work 
of adapting the mesh on the mesh generation program rather 
than on the actual numerical solution procedure of the gov- 
erning equations. In this method, a solution is first obtained, 
and regions of relatively large discretization errors are de- 
tected. A new mesh is then generated to concentrate points 
in regions where the large discretization errors occur. This 
new mesh may contain more or fewer points than the orig- 
inal mesh. 

For the second method, mesh movement, the number 
of points in the computational domain remains fixed. To im- 
prove the spatial accuracy of the solution, these points are 
moved into regions where solution gradients are relatively 
large. In general, this can be accomplished in two ways. The 
first way models the mesh as a spring network, where points 
are joined by linear springs with spring stiffnesses propor- 
tional to solution gradients. The mesh is then allowed to 
move into the relatively high gradient regions to produce ef- 
fectively a locally finer mesh. The second way uses forcing 
functions in a Poisson-equation grid generator to redistribute 
points. Either method of mesh movement is easily imple- 
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mented within existing solution algorithms because only the 
locations of the existing mesh points are changed. 

The final method of spatial adaptation is mesh enrich- 
ment. In this method points are added to regions of relatively 
large solution error by dividing locally the cells which make 
up the mesh or by embedding finer meshes in these regions. 
This method differs from mesh regeneration and movement 
in that the mesh is made finer in local regions while the 
global mesh topology remains the same. The method of 
mesh enrichment also is generally regarded as having ad- 
vantages over regeneration and movement, especially for 
transient problems. 2 For example, a distinguishing feature 
of mesh enrichment is that the original mesh is recovered 
once a refined feature has passed. Another feature is that 
the procedures are relatively fast compared to remeshing the 
entire domain. This feature is important since spatial adapta- 
tion is performed many times to track transient flow features. 

For the Euler and Navier-Stokes equations, computa- 
tional fluid dynamics algorithms are being developed based 
on spatial adaptation methods. With these equations, rel- 
atively large spatial discretization errors may be encoun- 
tered with flow features such as shock waves, shear layers, 
boundary layers, and expansion fans. These flow features 
can be resolved more accurately using the adaptation meth- 
ods mentioned above. Nakahashi and Deiwert, 3 for exam- 
ple, have used tension and torsional springs to move the 
mesh into regions where relatively large spatial discretiza- 
tion errors occur. This mesh movement approach showed 
considerable versatility for the problems treated. However, 
various constants were needed to control orthogonality and 
smoothness, and direct control of an optimal mesh adapta- 
tion procedure generally was not possible. Further exam- 
ples of spatial adaptation methods include the work of Usab 
and Murman. 4 In Ref. 4, embedded meshes of quadrilateral 
cells and nodes were used in regions of the mesh where 
shock waves occurred. This approach improved the spatial 
accuracy of the numerical method which resulted in highly 
accurate solutions for steady flow problems. Dannenhoffer 
and Baron 5 * 6 extended the work in this area using irregularly 
shaped embedded regions, which were coupled to the base 
mesh by a multiple-grid solution algorithm. Several other 
examples of spatial adaptation include methods which use 
flow solvers based on unstructured triangular and tetrahedral 
meshes in two and three dimensions, respectively. Peraire 
et al. 7,8 used mesh regeneration coupled with a finite ele- 
ment solution algorithm to sharply capture shock waves and 
complex shock structures. LOhner 9 developed a procedure 
in three-dimensions to locally enrich the mesh for transient 
flow problems by dividing tetrahedral elements which make 
up a base mesh to capture shock waves. Further, in this 
procedure, elements may be removed (coarsened) from the 
mesh if they are not necessary to produce a given level of 
spatial accuracy. More recently Kallinderis, Parthasarathy, 
and Wu 10 have developed and applied similar procedures. 


With respect to solution algorithms based on unstruc- 
tured meshes, the results published by the present authors 
demonstrated that these algorithms produce steady and un- 
steady solutions of comparable accuracy to results obtained 
using structured-grid solution algorithms. 11 In Ref. 1 1 struc- 
tured and unstructured mesh results were presented and com- 
pared for steady and unsteady flows where no mesh adap- 
tation was used. However, in Ref. 12, the present authors 
demonstrated that solutions of higher spatial accuracy are 
indeed possible for two-dimensional steady and unsteady 
flows through the use of mesh adaptation. Therefore, the 
purpose of this paper is to report on further modifications to 
the three-dimensional, unstructured-grid, upwind-type Euler 
code reported in Refs. 13, 14, and 15 to include mesh en- 
richment and coarsening procedures for steady and unsteady 
flow calculations. 

The objectives of the research are as follows: (1) to 
develop time-accurate mesh enrichment and coarsening pro- 
cedures for spatial adaptation, (2) to test the procedures by 
performing steady and unsteady calculations for a variety of 
cases, (3) to determine the accuracy of the spatially adapted 
solutions by making comparisons with published solutions 
produced by alternative methods and existing experimental 
data, and (4) to assess the efficiency of the spatially adapted 
solutions by making comparisons of required computer re- 
sources. 

The paper gives a detailed description of the mesh en- 
richment and coarsening procedures and gives a brief de- 
scription of the solution algorithm 13 ' 15 for completeness. 
Steady and unsteady flow results are presented to demon- 
strate an application of the adaptive mesh procedures. Steady 
flow results are presented for the ONERA M6 wing to as- 
sess the accuracy of the computed surface pressures by mak- 
ing comparisons with experimental data. Unsteady flow re- 
sults are presented in a three-dimensional simulation of a 
one-dimensional shock-tube problem to demonstrate an ap- 
plication of the enrichment and coarsening procedures for 
transient flows and to assess the accuracy of the computed 
results by making comparisons with the exact solution. 

Upwind-Type Euler Solution Algorithm 

The Euler equations are solved using the three- 
dimensional, unstructured-grid, upwind-type solution al- 
gorithm developed by Batina. 13 The solution algorithm 
of Ref. 13 was extended by the present authors 14, 15 for 
time-accurate, unsteady flow calculations on a deforming 
mesh. The algorithm, which is a cell-centered, finite- 
volume scheme, uses upwind differencing based on flux- 
vector splitting, 16 similar to upwind schemes developed for 
use on structured meshes. The flux-split discretization ac- 
counts for the local wave-propagation characteristics of the 
flow and captures shock waves sharply with at most one grid 
point within the shock structure. An additional advantage of 
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using flux-splitting is that the discretization is naturally dis- 
sipative and, consequently, does not require additional artifi- 
cial dissipation terms or the adjustment of free parameters to 
control the dissipation. However, in calculations involving a 
higher-order upwind scheme such as this, oscillations in the 
solution near shock waves are expected to occur. To elimi- 
nate these oscillations, flux limiting usually is required. In 
the present study, a continuously differentiable flux limiter 
was employed. 

The Euler equations are integrated in time using an 
implicit time-integration scheme involving a Gauss-Seidel 
relaxation procedure . 13 The relaxation procedure is imple- 
mented by re-ordering the elements that make up the un- 
structured mesh from upstream to downstream. The solu- 
tion is obtained by sweeping two times through the mesh 
as dictated by stability considerations. The first sweep is 
performed in the direction from upstream to downstream, 
and the second sweep is from downstream to upstream. For 
purely supersonic flows, the second sweep is unnecessary. 
This relaxation scheme is stable for large time steps and 
thus allows the selection of the step size based on the tem- 
poral accuracy of the problem being considered, rather than 
on the numerical stability of the algorithm. Consequently, 
very large time steps may be used for rapid convergence to 
steady state, and an appropriate step size may be selected for 
unsteady cases, independent of numerical stability issues. 

Spatial Adaptation Procedures 

In this section, the spatial adaptation procedures are de- 
scribed. These descriptions include detailed explanations of 
the procedures used to detect flow features and the proce- 
dures used to enrich and coarsen meshes of tetrahedra. 


Flow Feature Detection 

The first step of the spatial adaptation procedures is the 
detection of regions of relatively large discretization error 
so that the computational mesh can be locally enriched to 
improve the spatial accuracy or coarsened locally to reduce 
the computational costs. For numerical solutions of the 
Euler equations, these regions generally occur near flow 
features such as shock waves, stagnation points, slip lines, 
and expansion fans. The dominant flow feature for the cases 
considered in this study is a shock wave. 

There are a number of flow parameters that can be 
used for enrichment indicators based on the detection of 
shock waves. Parameters such as density, pressure, or total 
velocity are useful since these quantities are discontinuous 
through shocks. For example, first or second, divided or 
undivided differences in one of these parameters, similar to 
the work by Dannenhoffer and Baron , 6 can be used to detect 
shock waves. The enrichment indicator used in this study 


was the magnitude of the density gradient |Vp| which is 
often used to detect shock waves for steady flow problems. 

Mesh Enrichment 

Mesh enrichment of tetrahedral meshes is performed 
by starting with a relatively coarse mesh of cells and then 
subdividing these cells until a given level of spatial accuracy 
has been obtained. To prevent cells from being enriched too 
many times near flow discontinuities such as shock waves, 
an upper bound is placed on the number of times a cell 
can be divided. Presently, the upper bound is set before a 
calculation is performed where the upper bound is usually 
constrained by the computer memory available. For transient 
problems the mesh enrichment procedure may be performed 
at each time step of the integration of the governing flow 
equations or it may be performed once every set number of 
time steps. 

There are a number of ways to subdivide a tetrahedral 
cell. For example, a node could be added at the centroid 
of a cell and subdivided accordingly. This way, however, 
is not appealing because it often produces irregular cells, 
which tend to have an adverse effect on the accuracy of the 
solution algorithm. Another way is to add nodes at arbitrary 
locations along the edges of a cell. This is the approach 
taken in this study, but for convenience and to maintain 
regular cells, the nodes are added at the midpoints of the 
edges of the cells. 

Mesh enrichment is performed by using the enrichment 
indicator to determine if a cell is to be subdivided into 
smaller cells. To accomplish this, the enrichment indicator 
is computed for each cell and compared with a threshold 
value to determine whether a cell should be subdivided. In 
this study, the threshold value is set before the calculation 
is performed. If the threshold is exceeded, a new node is 
created at the midpoint of each edge of the tetrahedral cell, 
and the cell is subdivided into eight smaller cells. Special 
care must be taken, however, when an edge that is to be 
bisected lies on a boundary of the mesh, since the midpoint 
of the edge does not generally lie on the boundary. In this 
case, the location of the new node is determined by using 
a spline of the boundary coordinates. Further, the values 
of the flow variables for the new cells are determined by a 
linear interpolation of the conserved variables located at the 
nodes of the original cells. 

For a given tetrahedral cell to be enriched, either one 
edge, three edges (all part of the same triangular face), or 
all six edges are bisected. In the event that only two edges 
are marked to be bisected and they are on the same face of 
the cell, the third edge of the face is automatically bisected 
to prevent the creation of highly skewed or stretched cells. 
Similarly, if four or five edges of the cell are marked to 
be bisected, the remaining edges are bisected and the cell 
is fully enriched. Each time the mesh is enriched, a cell 
may be divided in one of three ways. The first way, shown 
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Fig. 1 Diagram of a type-8 element enrichment. 



(a) (b) (c) 

Fig. 2 Diagram of a type-2 element enrichment. 
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(C) 


Fig. 3 Diagram of a type-4 element enrichment. 


in Fig. 1, results when all six edges of a cell have been 
marked to be bisected. In this situation the cell is divided 
into eight new cells where the vertices of the inner cells are, 
in general, midpoints of edges that make up the original cell. 
The original cell is thus referred to as a type-8 element since 
after enrichment it becomes eight new tetrahedral cells. It 
should be noted that this way of subdividing a tetrahedron is 
not unique because there are three possible choices for the 
orientation of the inner diagonal edge that passes through the 
original cell. In the procedures, the orientation of the inner 
diagonal is based solely on the numbering of the nodes for 
the original cell being enriched. The second way, shown 
in Fig. 2, occurs if only one edge of the original ceil is 
marked to be bisected. In this situation, the marked edge 
is bisected, and two new cells are formed. The original 
ceil is thus referred to as a type-2 transition element since 
after enrichment it becomes two new cells. The third and 
final way, shown in Fig. 3, occurs if all three edges of a 
single face of the base cell are marked to be bisected. In 
this situation, the marked edges are bisected, and four new 
cells are formed. The original cell is thus referred to as a 
type-4 transition element since after enrichment it becomes 
four new cells. 

New cells formed from a type-8 element may be en- 
riched further. However, to prevent highly stretched cells, 
type-2 or type-4 transition elements are restricted from being 
divided further as indicated in Figs. 4 and 5. For cells from a 
type-2 transition element, if any of the nine edges that make 
up the two new cells are marked for enrichment, the original 
cell is made into a type-8 element as shown in Fig. 4. If in 
addition to this, either or both of the bottom two edges are 
marked (the lower left, right, or both), cells of the type-8 el- 



Fig. 4 Diagram illustrating details of further 


enrichment of type-2 transition elements. 

type-4 type-8 

element element 



Fig. 5 Diagram illustrating details of further 

enrichment of type-4 transition elements. 

ement are further enriched accordingly, as shown in Fig. 4. 
Similarly for cells from a type-4 transition element, if any 
of the fifteen edges that make up the four cells are marked 
to be bisected, the original cell is also made into a type-8 
element as shown in Fig. 5. If further, any edge that is part 
of the type-4 transition face is marked to be bisected, cells 
of the type-8 element are additionally enriched, as shown in 
Fig. 5. The type-4 transition face has nine edges that may 
be bisected further. If there is not a restriction on the edges 
marked to be bisected, a total of 512 permutations of further 
enrichment could result. However, since a triangular face 
of a type-4 transition element is prevented from having only 
two of its three edges marked for bisection, the total num- 
ber of possible permutations is reduced to 89. Eliminating 
permutations that are similar when rotated results in the 33 
enriched cells shown in Fig. 5. 
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(a) type-8 element coarsening. 



(b) type-4 transition element coarsening. 



(c) type-2 transition element coarsening. 

Fig. 6 Diagrams illustrating mesh coarsening possibilities. 

Mesh Coarsening 

Mesh coarsening of tetrahedral meshes is performed by 
removing added nodes and cells from previously enriched 
meshes to delete them from local regions of the mesh where 
certain flow features are no longer present. This procedure 
is necessary to adapt meshes to the numerical solution of the 
governing flow equations in order to minimize computational 
cost. Candidate cells for removal are cells that came from 
type-2, type-4, or type-8 elements and were marked for 
removal. Each time the mesh is coarsened, cells and nodes 
may be removed in one of several ways as shown in Fig. 6. 
For a type-8 element, three, five, or six nodes may be 
removed resulting in a type-4, a type-2, or an original cell, 
respectively (Fig. 6(a)). Similarly, if two of the three nodes 
that form the inner triangle of the face of a type-4 element 
are candidate nodes, the two nodes are removed and a type- 
2 element is formed as shown in Fig. 6(b). Likewise, if all 
three nodes that form the inner triangle of the face of a type- 
4 element are candidate nodes, the three nodes are removed 
and the one original cell that was divided previously into 
four remains (Fig. 6(b)). However, if only one of the three 
nodes that form the inner triangle of the subdivided face of 
a type-4 element is a candidate node then nothing is done. 
For a type-2 element there is only one node that may be 
removed which is the midpoint of a previously bisected edge. 
Removal of this node leaves only the cell that was divided 
originally into two as shown in Fig. 6(c). It should be noted 
that the mesh cannot become coarser than the original mesh. 


Results and Discussion 

Adaptive mesh results are presented in this section 
for a test case involving a simulated flow field using a 
mesh generated inside a cube, an ONERA M6 wing, and 
a three-dimensional simulation of a one-dimensional shock- 
tube problem. The results are used to assess the spatial 
adaptation procedures in three dimensions. The accuracy 
of the results are determined by making comparisons with 
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Fig. 7 Partial view of the surface of the original coarse 
mesh used in test case of the spatial adaptation 
procedures (1,805 nodes, 8,557 tetrahedra). 

results from alternative methods and available experimental 
data. 

Test Case for the 3D Adaptive Mesh Procedures 

To demonstrate the spatial adaptation procedures in 
three-dimensions a test case was performed. The test case 
was devised to assess newly developed procedures and data 
structures implemented within the mesh enrichment and 
coarsening procedures. The test case involves a mesh gen- 
erated inside a unit cube, where a partial view of the surface 
mesh is shown in Fig. 7. The mesh was generated using the 
advancing front mesh generation package, VGRID3D. 17 In 
Fig. 7 three sides of the surface mesh have been removed 
so that the interior of the mesh can be seen. The total mesh 
contains 1,805 nodes and 8,557 tetrahedra and serves as a 
starting mesh for the adaptive mesh procedures. 

The test case simulates an unsteady problem to demon- 
strate the mesh coarsening procedures. This case involves 
subdividing cells in the vicinity of the surface of a sphere, 
were the radius r of the sphere is increased with time at a 
constant rate. The equation for the surface of the sphere is 
given by 

x 2 + j/ 2 +z 2 = r 2 (1) 

where the origin of the sphere is the corner of the mesh. 
This equation is used to specify values of a field variable, 
\jj, for each cell in the mesh. For example, the value 
of 4> for a cell is zero if the centroid of the cell lies 
within ( x 2 + y 2 + z 2 < r 2 ) the sphere and the value of 
tp for a cell is one if the centroid of the cell lies outside 
(x 2 + y 2 + z 2 > r 2 ) the sphere. These values of 4> are 
used to compute the magnitude of the gradient of for 
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r = 0.50 (23,120 nodes, 129,795 tetrahedra). 


r = 0.75 (49,335 nodes, 281,820 tetrahedra). 


r = 1.00 (85,462 nodes, 489,369 tetrahedra). 

Fig. 8 Sequence of instantaneous surface meshes for the 
unsteady test case using three levels of enrichment. 



Radius 


Fig. 9 Variation of the number of tetrahedral cells in 
the mesh for the unsteady test case using 1, 

2, and 3 levels of enrichment. 

each cell in the mesh, where the magnitude of is 
used as the enrichment indicator. For the test case, three 
calculations were performed allowing one, two, and three 
levels of enrichment. For each calculation the spherical 
wave propagates in increments of A r = 0.025 of the radius 
for every time step and the calculation was continued until 
the radius r was 2.0. The sequence of surface meshes for the 
third calculation with three levels of enrichment is shown in 
Fig. 8 at three moments in time corresponding to a radius 
r of 0.50, 0.75, and 1.00. In this figure the surface meshes 
sharply define the spherical wave yet the adapted regions 
transition smoothly to the coarser regions of the mesh. A 
plot of the variation of the number of tetrahedral cells in 
the mesh versus the radius of the spherical wave is shown 
in Fig. 9 for the three calculations performed. In each 
calculation the number of ceils in the mesh varies smoothly 
as the spherical wave propagates through the mesh. This 
figure illustrates the rapid increase in the number of cells in 
the mesh as the enrichment levels are increased, which may 
be surprising for such a simple test case. 

ONERA M6 Wing 

The ONERA M6 wing was selected to further assess 
the adaptive mesh procedures in three dimensions. The ac- 
curacy of the adapted mesh results are assessed by making 
comparisons with experimental data and by making com- 
parisons with results obtained from other unstructured-grid 
Euler codes. The M6 wing has an aspect ratio of 3.8, a 
leading edge sweep angle of 30°, and a taper ratio of 0.56. 
The airfoil section is the ONERA D airfoil, which is a sym- 
metric 10% maximum thickness-to-chord ratio conventional 
section. The wing tip is rounded and is defined by a half 
body of revolution of the airfoil section. This wing has 
been widely studied and results have been obtained using 
many flow solvers on both structured 18 and unstructured 19 ' 21 
meshes. A steady flow calculation was performed at a 
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Fig. 10 Partial view of the original coarse surface mesh of 
the symmetry plane and the ONERA M6 wing. 
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Fig. 1 1 Steady-state convergence history for the ONERA 
M6 wing at M<» = 0.84 and a 0 = 3.06° 
using the adaptive mesh procedures. 

frees tr earn Mach number of 0.84 and an angle of attack 
of 3.06°. For this case, the starting mesh and the spatially 
adapted meshes are presented along with the corresponding 
density contour lines. Additionally, surface pressures were 
obtained for comparison with experimental data reported by 
Schmidt and Charpin. 22 

The mesh about the ONERA M6 wing configuration 
was generated using VGRID3D. The mesh extends 2^ wing 
semispans from the symmetry plane in the span direction. 
Also, the mesh extends 6^ root chordlengths above/below 
the wing surface as well as 6^ root chordlengths upstream 
and 10 root chordlengths downstream of the wing to rectan- 


upper surface lower surface 



(a) Original mesh. 


upper surface lower surface 



Fig. 12 Upper and lower surface meshes of the ONERA 
M6 wing for the original and adapted meshes. 


gular outer boundaries. The complete coarse mesh for the 
M6 wing contains 8,824 nodes and 46,5 16 tetrahedral cells. 
A partial view of the surface mesh for the symmetry plane 
and wing is shown in Fig. 10. 

The calculation was performed using implicit time inte- 
gration at a CFL number of 100,000. The final solution was 
obtained by adapting the mesh to the magnitude of the den- 
sity gradient every 300 iterations for the first 1,500 iterations 
and then marching the solution an additional 300 iterations 
on the final adapted mesh. Also, the adaptive mesh proce- 
dures allowed only one level of enrichment during the first 
900 iterations, and two levels were allowed thereafter. The 
convergence history for the calculation is shown in Fig. 1 1 
where the L 2 -norm of the density residual is plotted versus 
the CPU time in hours. In Fig. 1 1 the spikes in the L 2 -norm 
of the density residual mark the iterations at which spatial 
adaptation was performed. 

Figure 12 shows the upper and lower surface meshes 
for the original starting mesh and the final adapted mesh 
using two levels of enrichment. A summary of the sizes 
of the original and adapted meshes are given in Table 1. 
The information in the table consists of the total number of 
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upper surface lower surface 




Fig. 13 Upper and lower surface density contour 

lines (A p = 0.025) for the ONERA M 6 wing 
at Moo « 0.84 and a 0 ■ 3.06° using the 
original and adapted meshes. 

cells, nodes, boundary faces, wing boundary faces, and wing 
boundary nodes for each mesh. Notice in Table 1 that after 
one level of enrichment, the number of cells in the mesh is 
approximately twice as large as the number of cells in the 
original mesh. After two levels of enrichment, the mesh is 


Table 1 Comparisons of the ONERA M 6 wing meshes. 



Original 

Mesh 

One 

Level 

Two 

Levels 

Number of Cells 

46,516 

105,876 

352,417 

Number of Nodes 

8,824 

19,728 

63,855 

Number of Boundary 
Faces 

4,190 

7,606 

17,362 

Number of Wing 
Boundary Faces 

2,792 

6,048 

15,578 

Number of Wing 
Boundary Nodes 

1,419 

3,052 

7,822 



Fig. 14 Partial view of the adapted surface mesh of 
the symmetry plane and the ONERA M 6 
wing at Moo = 0.84 and a 0 = 3.06°. 

7.6 times larger. This increase in the number of cells is sig- 
nificant since it directly affects the computational resources 
required to perform a given calculation using the adapted 
meshes. Figure 13 shows the surface density contour lines 
on the original and adapted meshes using an increment of 
A p = 0.025. In Figs. 12 and 13 the upper surface is shown 
to the left and the lower surface is shown to the right. The 
upper surface contours (left part of Fig. 13) clearly show 
the lambda-shaped shock pattern formed by the two inboard 
shocks that coalesce near 87% semispan and then separate 
just outboard of 95% semispan (shown clearly for two levels 
of enrichment). The forward shocjc of the two shock region 
is a supersonic to supersonic shock. The aft shock of the two 
shock region is a supersonic to subsonic shock. The lower 
surface contours (right part of figure) indicate that there is 
little spanwise variation in density. A comparison of the two 
sets of contour lines reveals a considerable improvement in 
the resolution of the shocks when spatial adaptation is used. 
The meshes for the corresponding surface density contours 
show that points were clustered in the shock regions to pro- 
duce shock waves that are spatially sharp. Also, points were 
clustered near the leading edge to improve the accuracy of 
this high gradient region. A partial view of the surface mesh 
and density contour lines for the symmetry plane and wing 
are shown in Figs. 14 and 15. In Fig. 14 the mesh in the 
symmetry plane gives an indication of the mesh spacing off 
of the surface of the wing. Similarly, in Fig 15 the con- 
tour lines give an indication of the spatial resolution of the 
solution off of the surface of the wing. 

To assess the accuracy of the results, the calculated sur- 
face pressure coefficients are compared to the experimental 
data at six semispan stations for the original and adapted 
meshes. The semispan stations are at 77 = 0.20, 0.44, 0.65, 
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Fig. 15 Partial view of the surface density contour lines 
(A p = 0.025) on the symmetry plane and on the 
ONERA M6 wing at Moo = 0.84 and a 0 = 3.06° 


0.80, 0.90, and 0.95. Here the upwind-biased flow variables 
of surface triangles, with a common edge along the semis- 
pan station, are interpolated to the edge to determine the 
surface pressures. The comparisons of surface pressures are 
shown in Figs. 16 and 17 for the original and the adapted 
mesh involving two levels of enrichment. In these figures 
the calculated Euler results are given by solid curves where 
plus symbols have been used to mark the interpolated values 
along the semispan stations. The experimental data is de- 
noted by circle and square symbols representing upper and 
lower surface pressure coefficient data. A comparison of 
the calculated surface pressures in Figs. 16 and 17 shows an 
improvement using the adapted mesh (Fig. 17), especially 
on the upper surface near the leading edge where the suc- 
tion peak was poorly predicted in Fig. 16 using the coarser 
mesh. For Fig. 17 at ij = 0.20 there are two shock waves on 
the upper surface. The forward shock is well predicted and 
agrees well with the data. The aft shock wave is too strong 
and located too far back on the wing when compared with 
the data. This however is consistent with other inviscid flow 



0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 


x/c 


x/c 


Fig. 16 Comparison of calculated and experimental values 
of pressure coefficient computed using the 
original coarse mesh for the ONERA M6 
wing at Moo = 0.84 and ao = 3.06°. 




0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 


x/c x/c 

Fig. 17 Comparison of calculated and experimental values 
of pressure coefficient computed using the 
adapted mesh (2 levels) for the ONERA M6 
wing at Moo = 0.84 and ao = 3.06°. 


9 




calculations for this case. At rj = 0.44 the forward shock has 
moved aft and the aft shock has moved forward. Again the 
forward shock wave is well predicted and the aft shock is 
slightly strong in comparison with the data. At 77 = 0.65 the 
forward shock is near 19% chord and the aft shock wave is 
near 50% chord. Both shock waves are well predicted, and 
the pressure level between the two shock waves agrees well 
with the data. At 7? = 0.80 the shock waves have begun to 
coalesce near 30% chord. A feature of this data worth not- 
ing is that the spatial adaptation procedures have helped to 
clearly define the two shock waves at 77 = 0.80. Calculations 
by other researchers 18-20 have failed to isolate the two shock 
waves at this semi-span location mainly due to the coarsness 
of the meshes used. At 77 = 0.90 the two shock waves have 
merged to form a single strong shock near 30% chord. The 
shock wave is sharply captured with one interior grid point, 
which is common for upwind schemes of this type. Finally, 
at r} = 0.95 the upper surface pressure results show a strong 
shock slightly aft of the experimental data. Also, just aft 
of the leading edge, the upper surface results slightly un- 
derpredict the data. In general the lower surface pressure 
coefficients agree well with the data at all semispan stations. 

The computed results are further assessed by comparing 
the lift coefficients obtained using the original mesh and the 
adapted mesh (two levels of enrichment) with those reported 
in Refs. 20, 23, and 24. This comparison is shown in 
Table 2, where the lift coefficient in Ref. 20 was computed 
on a mesh with 53,989 nodes and 288,170 cells, the lift 
coefficient in Ref. 23 was computed on a mesh with 16,984 
nodes and 231,507 cells, and the lift coefficient in Ref. 24 
was computed on a mesh with 173,412 nodes and 1,013,718 
cells. For these values of lift coefficient a reference area of 
S rt j = 0.5255 was used. The table shows good agreement 
between the lift on the final adapted mesh (two levels) and 
the lift obtained using the other unstructured Euler solvers 
giving confidence in the adapted mesh solution. 

The calculation for the final adapted mesh solution was 
run for 1,800 iterations and required approximately 12.5 
hours of CPU time and 125mw of memory on the Cray- 
2 computer at the NAS facility located at the NASA Ames 
Research Center. 


Table 2 Comparisons of lift coefficients for the ONERA 
M6 wing at M ^ = 0.84 and a 0 = 3.06°. 



Original 

Two 

Ref. 

Ref. 

Ref. 


Mesh 

Levels 

20 

23 

24 

Lift Coef. 

0.2827 

0.2901 

0.2923 

0.2911 

0.2901 


Shock-Tube Problem 

An unsteady one-dimensional shock-tube problem was 
used to evaluate the accuracy and efficiency of the spatial 
adaptation procedures in three-dimensions. The shock-tube 
problem is illustrated in Fig. 18 where a diaphragm sepa- 



<« I | | | [ O) ] (3) j (4) 


Fig. 18 Illustration of a shock-tube problem. 

rates a high pressure (compression) chamber (1) and a low 
pressure (expansion) chamber (4). Initially the pressure dis- 
tribution in the shock-tube is an ideal “step”. At the instant 
the diaphragm bursts the initial pressure “step” separates 
into a shock wave, which propagates to the right into the 
low pressure chamber, and an expansion fan, which propa- 
gates to the left into the high pressure chamber. The region 
traversed by the shock (3) and the region traversed by the 
expansion fan (2) is separated by a contact surface. Each in- 
terface between the four regions moves at a constant speed 
as shown in Fig. 18. 

The shock-tube problem provides a good test case for 
the adaptive mesh procedures since the flow contains many 
features expected to occur in transient problems such as a 
moving shock wave, an expansion fan, and a contact surface. 
Also, the exact solution 25 is available for comparison, and 
the computational resources for this problem are relatively 
small. The problem is a challenging one for the spatial 
adaptation procedures since all the flow features must be 
tracked accurately in time as the solution progresses. 

The test case considered involves air at the same tem- 
perature in the low and high pressure chambers and a di- 
aphragm pressure ratio of five. The initial values of the 
primitive variables on each side of the diaphragm are given 
below in Eq. 2. 

*=5 

p 4 

m = v\ = u>i = 0 

( 2 ) 

ti 4 = 1>4 = W4 = 0 

P 4 

The calculations were performed on a coarse mesh, within 
parallelepiped boundaries, of unit length and with width and 
height that is four percent of the length of the tube. The 
mesh was generated using the VGRID3D mesh generator. 
Figure 19 shows the surface mesh for the parallelepiped. 
The total mesh contains 1,800 tetrahedra and 562 nodes. 

The first calculation for the shock-tube problem was 
performed on the coarse mesh to obtain a solution for com- 
parison with the exact solution. The coarse mesh result 
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Fig. 19 Surface mesh for the parallelepiped. 

also provides a solution for comparison with the spatially 
adapted results. The solution was obtained using implicit 
time-marching with a nondimensional time step of 0.001. 
The calculated density profiles at three moments in time, 0.1, 
0.2, and 0.3 are shown in Fig. 20. These profiles were ob- 
tained by sorting all 1,800 cell-centered values of the prim- 
itive variables according to their x-coordinate locations and 
plotting every other point. Since the three-dimensional so- 
lution is plotted as a function of the x-coordinate there may 
be a variation of the solution in the y- and z-coordinate di- 
rections. As the solution in Fig. 20 progresses in time, the 
expansion fan spreads and moves to the left while the shock, 
followed by the contact surface, moves to the right. Figure 
21 shows the surface density contour lines at the same three 
moments in time. The second calculation used the spatial 
adaptation procedures starting with the coarse mesh of the 
first calculation. The mesh was locally pre-embedded about 
the initial pressure discontinuity for two levels of enrichment 
in order for the initial discontinuity to be sharply defined. 
The solution was marched in time for 3,000 time steps and 
the mesh was adapted to the magnitude of the density gradi- 




Fig. 20 Comparison of calculated and exact density 
profiles in a shock-tube for a sequence of 
times (0.1, 0.2, 0.3) on the original coarse 
mesh (every other point plotted). 



Fig. 21 Surface density contour lines (A p = 0.0067) for a 
shock-tube for a sequence of times (0.1, 0.2, 

0.3) using the original coarse mesh. 

ent once every 20 time steps using a threshold value of 1 .0. 
During the course of the calculation, mesh refinement was 
restricted to two levels of enrichment. Figure 22 shows 
the resulting density profiles using spatial adaptation for the 
same three moments in time that were shown previously. 
Similar to the previous results, these profiles were obtained 
by sorting all the cell-centered values of the primitive vari- 
ables according to their x-coordinate locations and plotting 
every 20th point. In this figure the amount of adaptation can 
be seen in the distribution of the calculated data throughout 
the tube, where the majority of the cells are concentrated in 
the regions of the expansion fan, the contact surface, and the 
shock wave. The concentration of cells can be seen more 




0.0 0.2 0.4 0.6 0.6 1.0 


Fig. 22 Comparison of calculated and exact density 
profiles in a shock-tube for a sequence of 
times (0.1, 0.2, 0.3) on the adapted mesh 
(every 20th point plotted). 
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Fig. 23 Surface mesh and density contour lines 

(A p = 0.0067) for a shock-tube for a sequence of 
times (0.1, 0.2, 0.3) using an adapted mesh. 


clearly in Fig. 23 where the surface mesh and surface density 
contour lines are shown. The figure illustrates, by compar- 
ison with Fig. 21, that the shock wave is sharply captured 
using the spatial adaptation procedures. Figure 22 shows 
good agreement in temporal and spatial accuracy with the 
exact solution giving confidence in the adaptive mesh pro- 
cedures. The shock-tube results obtained using the spatial 
adaptation procedures required approximately 5.5 hours of 
CPU time and 18.2mw of memory on the Cray-2s computer 
at NASA LaRC. Also, the spatial adaptation procedures con- 
sumed approximately 20% of the total CPU time. 

Concluding Remarks 

Spatial adaptation procedures for the accurate and effi- 
cient solution of steady and unsteady inviscid flow problems 
were described. The adaptation procedures were developed 
and implemented within a three-dimensional unstructured- 
grid upwind-type Euler code. These procedures involve 
mesh enrichment and mesh coarsening to either add points 
in high gradient regions of the flow or remove points where 


they are not needed, respectively, to produce solutions of 
high spatial accuracy at minimal computational cost. 

Steady and unsteady results were presented to demon- 
strate applications of the spatial adaptation procedures to 
three-dimensional problems. Steady transonic flow results 
were obtained for the ONERA M6 wing to assess the ac- 
curacy of the computed surface pressures by making com- 
parisons with experimental data. The results obtained us- 
ing spatial adaptation were found to be in good agreement 
with the experimental data giving confidence in the mesh 
enrichment and coarsening procedures. Unsteady flow re- 
sults were obtained in a three-dimensional simulation of a 
one-dimensional shock-tube problem to demonstrate an ap- 
plication of the mesh enrichment and coarsening procedures 
for a time-dependent problem. The accuracy of the com- 
puted results was assessed by making comparisons with the 
exact solution. Both the steady and unsteady solutions ob- 
tained using spatial adaptation were shown to be of high 
spatial accuracy, primarily in that the shock waves were 
sharply captured. 
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